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Previous experiments have revealed that shock waves driven through dissipative gases may become 
unstable, for example, in granular gases, and in molecular gases undergoing strong relaxation effects. 

The mechanisms controlling these instabilities are not well understood. We successfully isolated and 
investigated this instability in the canonical problem of piston driven shock waves propagating into 
a medium characterized by inelastic collision processes. We treat the standard model of granu¬ 
lar gases, where particle collisions are taken as inelastic with constant coefficient of restitution. 

The inelasticity is activated for sufficiently strong collisions. Molecular dynamic simulations were 
performed for 30,000 particles. We find that all shock waves investigated become unstable, with 
density non-uniformities forming in the relaxation region. The wavelength of these fingers is found 
comparable to the characteristic relaxation thickness. Shock Hugoniot curves for both elastic and 
inelastic collisions were obtained analytically and numerically. Analysis of these curves indicate 
that the instability is not of the Bethe-Zeldovich-Thompson or Dyakov-Kontorovich types. Analysis 
of the shock relaxation rates and rates for clustering in a convected fluid element with the same 
thermodynamic history outruled the clustering instability of a homogeneous granular gas. Instead, 
wave reconstruction of the early transient evolution indicates that the onset of instability occurs 
during the re-pressurization of the gas following the initial relaxation of the medium behind the lead 
shock. This re-pressurization gives rise to internal pressure waves in the presence of strong density 
gradients. This indicates that the mechanism of instability is more likely of the vorticity-generating 
Richtmyer-Meshkov type, relying on the action of the inner pressure waves development during the 
transient relaxation. 
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I. INTRODUCTION 

Shock waves driven into dissipative gases sometimes 
develop instabilities. Granular media, which are char¬ 
acterized by inelastic particle collisions, is one example. 
Previous experiments have identified unstable formations 
of finger-like jets in granular media dispersed by shock 
waves driven through air DU and for rapid granular 
flows down a chute [3]. Similar pattern formations can 
be seen when granular media are subjected to a verti¬ 
cally oscillating bed, both experimentally and numeri¬ 
cally DD- In the latter, the periodic agitation of the 
container walls drive strong shocks and expansion waves 
into the non-uniform granular gas. The complex tran¬ 
sient dynamics involved in the these past configurations 
have prevented the authors to clearly identify the mech¬ 
anisms controlling the instability. In the present study, 
we pose the problem in the classical formulation of a pis¬ 
ton suddenly accelerated to a constant speed into a gas 
medium, as illustrated in Figure [T] 

Previous investigations of this canonical problem have 
looked at the one-dimensional structure and evolution of 
shock waves in granular gases, although instabilities had 
not been identified [HID- Goldshtein et al. revealed that 
the structure of shock waves driven by a piston into a 
granular gas is composed of three distinct regions [Sj . The 
first region follows the shock front, and is composed of a 
rapid increase in granular temperature (region I). Due to 
the inelasticity and increased rate of the collisions within 


this excited region, the granular temperature of the ma¬ 
terial falling further behind the shock starts to decrease, 
while density increases; this marks the ‘relaxing’ region 
(region II). Eventually, the collision amplitudes become 
sufficiently weak such that visco-elastic particles collide 
elastically. In this ‘equilibrium’ region (region III), the 
gas retains a finite granular temperature. When all colli¬ 
sions are assumed inelastic, the equilibrium region tends 
to zero granular temperature. 

Kamenetsky et al. [7] investigated the evolution of such 
a structure numerically by solving the one-dimensional 
Euler equations for granular media. The authors revealed 
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FIG. 1: Temperature distribution for a thermally 
relaxing shock wave travelling at velocity D 
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interesting dynamics prior to the shock wave attaining 
the developed structure illustrated in Figure |T| In par¬ 
ticular, the authors found that the lead shock front pulls 
back towards the piston for a short period, before attain¬ 
ing a constant velocity. The dynamics of this stage were 
not explained nor further explored. Nevertheless, as we 
will show in the present article, these turn out to have a 
strong bearing on the multi-dimensional shock instabil¬ 
ity. 

Qualitatively, a structure similar to that shown in 
Figure [T] is observed for sufficiently strong shock waves 
driven into molecular gases, whereby the shock is 
strong enough to bring about inelastic collisions between 
molecules (i.e., via endothermic reactions) [5j. Interest¬ 
ingly, these types of relaxing shock waves have also been 
shown to sometimes become unstable. Unstable shock 
structures have been observed experimentally in suffi¬ 
ciently strong shocks leading to ionization Emu, disso¬ 
ciation m and in gases with high specific heats EnHHI- 

Current models for predicting such shock instability 
are mostly based on jump conditions between the initial 
and final equilibrium states, without knowledge of the 
kinetic processes linking the two states. The D’Yakov- 
Kontorovich (DK) and the Bethe-Zeldovich-Thompson 
(BZT) mechanisms require the shape of the Hugoniot 
curves to have anomalous properties (see, for example, 
Refs. [5] and m)- The Hugoniot curve is the lo¬ 
cus of the equilibrium post shock state, usually repre¬ 
sented in the pressure-specific volume plane. While the 
Hugoniot curves can be obtained experimentally for a 
given substance, investigation of their properties in the 
context of BZT and DK instabilities predicted stable 
shocks at experimental conditions corresponding to un¬ 
stable shocks mm- 

Another mechanism of interest involved in shock in¬ 
stability is that of Richtmyer-Meshkov and Rayleigh- 
Taylor type instabilities, although such instabilities have 
not been reported in the cases above. In such a multi¬ 
dimensional instability, misaligned gradients of density 
and pressure lead to vorticity production |X6 ]. This type 
of instability is a universal physical phenomena encoun¬ 
tered, for example, in gases B7I, plasmas ^TSj , Bose- 
Einstein condensates m , and combustion m- 

Models for predicting the instability of relaxing shocks 
involving the kinetics of the relaxation process have only 
very recently been formulated. Direct numerical simula¬ 
tions at the continuum level in the case of ionizing shocks 
has indeed recovered the instability EH 122, suggesting 
that it is related to the hydrodynamic coupling with the 
kinetics of the relaxation process. This suggests that an 
account for the kinetics of the relaxation process may 
be required to predict the shock instability in relaxing 
media. 

In the absence of bulk flow, it has been shown that 
such dissipative gases are subject to clustering instabil¬ 
ities (23H25I . This clustering instability , first shown by 


Goldhirsch and Zanetti [23], is seen in granular gases, 
where the collisions can be assumed to remain inelastic 
for all impact conditions. In such a medium, an initially 
homogeneous gas develops clusters during its cooling, 
which takes the form of filamentous structures. Gas is 
preferentially accelerated towards regions of higher den¬ 
sity, owing to the local greater rate of pressure decay in 
these regions due to dissipation. Since the material pass¬ 
ing through the shock structure undergoes the same cool¬ 
ing process, the clustering instability may be controlling 
the local non-homogeneities within the shock structure. 
This link is further explored in the present paper. 

To summarize, the goal of our present study is twofold. 
We wish to first isolate the shock instability in relaxing 
media in a canonical problem, conducive to further anal¬ 
ysis. Second, we wish to determine the mechanism con¬ 
trolling the instability. The qualitative correspondence 
of the structure of granular gases and molecular gases 
suggests that both problems can be studied by the same 
formalism, provided the collision properties are modified 
to account for the finite temperature equilibrium region 
of molecular gases. 

To investigate the evolution and stability of such shock 
waves, we adopt the simple kinetic model previously used 
to describe dissipative granular gases by Goldhirsch and 
Zanetti: the collision between “hard particles” of finite 
radius is modeled deterministically using a constant co¬ 
efficient of restitution taken below unity. This model is 
the simplest kinetic model that can mimic relaxation. 
In order to capture the structure of relaxing gases more 
closely, we also assume that the collisions are activated 
by an impact energy threshold. Such a threshold is also 
applicable to granular media, which has been used to 
better imitate the visco-elastic behavior of colliding par¬ 
ticles [25] . 

The paper is organized as follows: Section II outlines 
the details of the molecular dynamics model used in this 
study. Section III addresses the evolution and struc¬ 
ture of shock waves predicted by the molecular dynamic 
model. Section IV provides further discussion and anal¬ 
ysis of the mechanism controlling the shock instability. 
Finally, Section V offers our closing remarks. 

II. DETAILS OF MOLECULAR DYNAMICS 
MODEL 

The approach we use is a deterministic hard particle 
dynamic approach in a 2D environment, akin to the prob¬ 
abilistic approach of the Direct Simulation Monte Carlo 
(DSMC) technique [27] ■ In such models, only the col¬ 
lision rules are prescribed in order to capture a phys¬ 
ical phenomenon (granular gases, relaxation, chemical 
reactions, etc.). We employ the standard deterministic 
method used for granular gases, both in its kinetic theory 
and in particle based simulations. Indeed, much of the 
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kinetic theory of dilute, idealized gases can be obtained 
by treating molecules as hard spheres with no internal 
structure [28] [29] . 

The current model assumes that collisions with bound¬ 
aries are elastic, yielding a symmetry condition that is 
implemented in order to not artificially introduce sup¬ 
plementary system size effects. Each binary collision is 
elastic, unless an activation threshold is reached. The 
post-collision velocities of two particles are calculated as: 


<5 =*—(! +035 

u'j =Uj + -(1 + s*)g?j 


(1) 


where = u™ — Uj is the normal component of the 
relative velocity of the two disks. 

Activation is assumed to occur when the collision be¬ 
tween two disks is sufficiently strong. This mimics the 
excitation of higher degrees of freedom (rotation, vibra¬ 
tion, dissociation, ionization, etc.) with increasing tem¬ 
peratures [25]. This is also a simple model for granular 
media undergoing visco-elastic collisions [25]. Quantita¬ 
tively, the collision between two disks is assumed to be 
elastic if g-‘j is below a threshold u*, a classical activa¬ 
tion formalism in chemical kinetics. For collisions with 
a higher amplitude, we assume an inelastic dissipative 
collision, which is modeled with a constant coefficient of 
restitution e < 1 . i.e., 


r i ^ m<u* 

\e if \g^\>u* 


( 2 ) 


where the predefined u* and £ remain constant during 
each simulation. 

The problem we study is a classical shock propagation 
problem, whereby the motion of a suddenly accelerated 
piston driven in a thermalized medium drives a strong 
shock wave. The driving piston is initially at rest and 
suddenly acquires a constant velocity u p . Collisions with 
the piston are elastic. This model allows for the dissipa¬ 
tion of the non-equilibrium energy accumulated within 
the shock structure, which terminates once the collision 
amplitudes fall back below the activation threshold. In 
this manner, the activation threshold also acts as a tun¬ 
able parameter to control the equilibrium temperature in 
the post shock media. Note that the model assumed is 
also the standard model for granular gases m , allow¬ 
ing us to compare with the established hydrodynamic 
description of this type of media. 

The MD simulations thus reconstruct the dynamics of 
hard disks. These are calculated using the Event Driven 
Molecular Dynamics (EDMD) technique first introduced 
by Alder and Wainright [31] ■ We use the implementation 
of Poschel and Schwager ;32j, that we have extended to 
treat a moving wall (piston). The particles were initial¬ 
ized with equal speed and random directions. The sys¬ 
tem was let to thermalize and attain Maxwell-Boltzmann 


statistics. Once thermalized, the piston started moving 
with constant speed. This code was implemented and 
tested for non-dissipative media in our previous study 
[35], where the simulated shock jump conditions agreed 
with those which were derived for hard disk mixtures. 

The initial packing factor of the disks was chosen to be 
r ]i = (Nnd 2 )/4:A = 0 . 012 , where NttcP/A is the volume 
(area) of the N hard disk with diameter d, and A = 
L x x L y the domain area; the initial gas is thus in the 
ideal gas regime [33]. All distances have been normalized 
by the initial mean free path of the system of disks Ai, 
which takes the form [30] : 


Ai 


1 

V / 27rdni6 2 (??i) 


(3) 


where 62 (^ 7 ) = (1 — y^)/(I — rj) 2 is the Enskog factor 
for a 2D system of hard particles, and n\ = N/A is the 
initial number density of particles. All speeds are scaled 
by the initial root mean squared velocity w rms (i) of the 
disks, fixing the time scaling by the initial mean free time 
D = Ai/u rms (i). 

The numerical experiments were performed using 
30,000 disks, unless otherwise noted. A domain size of 
L x x L y = 172.9 x 17.2 and disk radius <7 = 0.019 was 
used to satisfy the packing factor of r]i = 0.012. The di¬ 
mensions of the domain, with 30,000 particles, was found 
to be an appropriate size to investigate and capture in¬ 
stability, allowing for sufficiently fast computing in order 
for results to be ensemble averaged. Ensemble and coarse 
grain averaging was implemented to investigate the one- 
dimensional shock structure. For each set of parameters, 
an ensemble of 50 simulations was taken, with the macro¬ 
scopic properties taken in strips of width Ax ~ 0.5Ai 
parallel to the piston face. 

All macroscopic properties are scaled by the initial 
state, unless otherwise noted. The density p is taken 
by tracking the number of disks within each strip, and 
the granular temperature is taken with the root mean 
squared velocity, i.e., T = ^u 2 ms . The pressure is ap¬ 
proximated from the Helfand equation of state for elastic 
disks [33 ] : 


pT 

P (1-t?) 2 


(4) 


To investigate the dynamics of the shock waves, the 
family of characteristics were constructed. The particles 
paths (P), forward (C + ) and backward ( C ~) running 
characteristics on an x vs. t plane are given by: 

_ dx v dx + _ dx _ 

P : =u C + : — r±- =u + c C : —— =u-c 

dt dt dt 

(5) 

where u is the local particle velocity normal to the piston 
and c is the local speed of sound, at a given time. They 
represent the trajectories of fluid particles, right running 
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pressure waves and left running pressure waves, respec¬ 
tively |T5] . The scaled speed of sound for such a media is 
approximated for an elastic system of disks, taken as {33] : 

+ + 277(1-77)-!) (6) 

^rms i V J-1 

The local packing factor is taken from the density jump, 
V = mp/pi- 

The trajectories of the characteristics were obtained 
numerically by integrating ([5]). The C + characteristics 
are initiated from the piston face at specified intervals 
in time, while C~ characteristics are initiated from the 
shock front at similar time intervals. Particle paths are 
initialized at specified locations away from the initial pis¬ 
ton position, denoted as £ = x(t = 0) for each path. 

III. RESULTS 

In this section we discuss the results obtained using 
the described model. We compare the evolution of shock 
structure and ensuing instability for varying properties. 
First, we look at the evolution of shock structure in detail 
for a single case. Next, we perform a parametric study 
to see how the evolution, shock structure, and stability 
vary with u p , u* and e. 

A. Evolution of shock structure 

The first case we look at is for u p = 20, u* = 10, and 
e = 0.95. Figure [2] shows an example of the evolution 
of the one-dimensional temperature distribution. In ad¬ 
dition to showing the instantaneous structure, the peak 
temperature and temperature at the piston are tracked. 
Initially the temperature jump of the shock is approx¬ 
imately ~ 400, as predicted for a system of elastic 
disks [331]. The temperature measured at the piston sur¬ 
face decays until coming to a quasi-equilibrium state, at 



FIG. 2: Evolution of shock structure for u p = 20, 
u* = 10, and £ = 0.95. 


which point most inelastic collisions have subsided - note 
that the kinetic model taken maintains an exponentially 
small fraction of activated collisions as the temperature 
decays below the activation temperature. The peak tem¬ 
perature also decays initially, which is followed by an 
oscillation before reaching an equilibrium peak tempera¬ 
ture. These dynamics are very similar to the ones pre¬ 
dicted by Kamenetsky et al. in inviscid hydrodynamic 
simulations of granular gases with a constant e [7j. 

Figure [3] shows the evolution of the averaged temper¬ 
ature, density and pressure fields in the x - t plane, in a 
frame of reference moving with the piston. Selected par¬ 
ticle paths, C + characteristics extending from the piston, 
and C~ from the shock front are also shown in order to 
more clearly illustrate the dynamics. For example, the 
shock is the locus along which all forward facing pressure 
waves C + coalesce. The shock wave driven by the piston 
generates an increase in the medium pressure, density 
and temperature. As the medium behind the shock be¬ 
gins to cool, the lead shock is seen to decay. The cooling 
of the gas and decay of the lead shock can be correlated 
by the forward facing pressure waves. The excess relax¬ 
ation behind the lead shock leads to an eventual pull¬ 
back of the lead shock towards the piston. A similar 
pullback was observed by Kamenetsky et al. in their hy¬ 
drodynamic simulations [7]. 

The cooling of the gas behind the lead shock, which 
can be followed along the corresponding particle paths, 
eventually is punctuated by an increase of density and a 
re-pressurization. This can be clearly observed at t ~ 2. 
The origin of this re-pressurization is not clear at present, 
but may be correlated with the arrival of the rear fac¬ 
ing pressure waves (along the C~ characteristics shown), 
originating at the decaying shock. Interestingly, the rear 
re-pressurization leads to a forward-facing pressure wave, 
arriving at the lead shock at t ss 3. This marks the re¬ 
acceleration of the lead shock towards its final equilib¬ 
rium structure. 

Figure [4] shows the evolution of shock morphology for 
this case, obtained from a single realization. These re¬ 
sults show the birth of an unstable structure, which we 
distinguish by density perturbations and corrugations ap¬ 
pearing within the shock structure. Initial stages of the 
evolution do not show distinguishable instabilities, as 
seen at t = 0.3, and up to t = 1.5. This is the point 
where the shock front stops propagating ahead of the 
piston. For later times, instabilities in the form of high 
density clusters and corrugations appear at the piston 
face. This is seen at t = 2.7, confirming that these insta¬ 
bilities occur between t = 1.5 and t = 2.7. 

Comparing with the evolution of pressure shown in 
Figure [3]{c), this range in time is when the early particle 
paths undergo a re-pressurization event on route to at¬ 
taining an equilibrium state. This indicates that the ori¬ 
gin of the instability may be associated with this distinct 
feature of the relaxation process; a possible mechanism is 
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FIG. 3: (Colour online) Evolution of (a) temperature, (b) density and (c) pressure on an x vs. t plane, in the piston 
frame of reference, for u p = 20, u* = 10 and e = 0.95. Evolutions shown with select particle paths (solid white), 
forward (solid black) and backward (dashed blue) running characteristics. 



FIG. 4: Evolution of shock morphology, for a single realization, at t = (a) 0.3, (b) 1.5, (c) 2.7 and (d) 3.9; with 

u p = 20, u* = 10 and e = 0.95. 


discussed in Section IV. Once the shock evolution enters 
the developed stage, the clusters begin growing from the 
piston, as demonstrated by the snapshot at t = 3.9. 


Figure [5] shows the particle distribution in the shocked 
material in relation to the mean temperature and density 
distributions. Superposed on the particle distribution 
plot is the coarse grained velocity vector field. This in¬ 
stantaneous vector field is rendered using streamlines, in 
order to better visualize the existence of coherent struc¬ 
tures. The streamlines were obtained by interpolating 
on the uniform grid of coarse grained averaged velocity 
vector field. Results show that substantial disturbances 
in speed are present in the region of the high density 
gradients. Streamlines converge toward the high density 
fingers, giving rise to convective rolls. 


B. Parametric study of the shock structure and its 
evolution 

Dimensional analysis and independent parameters 

The macroscopic dynamics of the model introduced is 
expected to have a relatively small number of controlling 
parameters. Dimensional analysis permits us to deter¬ 
mine the number of parameters controlling the dynamics. 
The initial thermodynamic state is uniquely defined by 
its granular temperature Ti, density p\ and packing frac¬ 
tion r)i. The shock dynamics, depend on the piston speed 
Up , the activation threshold u* and the degree of inelas¬ 
ticity, e. Furthermore, we are interested at conditions in 
which the strong shock limit applies and the initial in¬ 
ternal energy does not control the dynamics 0GSj; this 
is the case where the experimental observations of shock 
instability have been made, for both the granular and re¬ 
laxing molecular gases, as discussed in the Introduction. 
Under the scaling of our variables, this reduces to the 
limit where u p 1 and u* 1. Under this limit, the 
parameters of the problem reduce to u p /u *, rj, and s. 
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FIG. 5: Particle distribution and coarse-grained 
stream-lines for a single realization (top), with ensemble 
and coarse grained distributions of temperature and 
density (bottom) after t = 8.13, with u p = 20, u* = 10, 
and £ = 0.95. 


Dependence on u v /u * 


It was found that u p /u* controlled the type of dynam¬ 
ics observed during the relaxation process. Figure[7]com- 
pares the evolution of the temperature and pressure fields 
obtained for u p /u* = 1.0, 1.5, and 2.0, with £ = 0.95. For 
the case u p /u* = 1.0 shown in Figure [7](a). In this case, 
the strong initial shock wave is followed by a gradual de¬ 
cay of the shock velocity. This decay does not cause the 
shock to pull back towards the piston, and the early par¬ 
ticle paths do not experience a re-pressurization along 
the piston face. When u p /u* (Figure [7](b)) is increased 
to 1.5, the shock front stalls with respect to the piston 
and a moderate re-pressurization is seen along the piston 
face. Further increase of the piston speed leads to a more 
marked shock pull-back and re-pressurization event, such 
as that seen in Figure[7](c) for u p /u* = 2.0. The threshold 
for oscillatory behavior for the front shock and internal 
re-pressurization is approximately u p /u* 


1 . 


Figs. 8(a) and |8(b)| show a comparison of the developed 
distributions of density and kinetic energy after equal 
piston displacement for u p /u* = 1.0, 1.5, and 2.0, with 
£ = 0.95. Both distributions show that the distance of 
the shock front decreases as u p /u* increases. This is at¬ 
tributed to the decreasing relaxation zone length for in¬ 
creasing u p /u*, which is seen by the steeper slopes for in¬ 
creasing density and decreasing kinetic energy. The peak 
energy increases with increasing u p /u*, as expected for 
increasing u p . All cases share a common kinetic energy 
at the piston face, corresponding to the quasi-equilibrium 
state with kinetic energy tending to 5-8% of the activa¬ 
tion energy. 


In order to validate the results of our dimensional anal¬ 
ysis, we varied both u p and u* in the hypersonic limit 
u p ^ u* 1. Figure 6(a) and |6(b) shows results for 
the distributions of density and kinetic energy, respec¬ 
tively, after equal piston displacement while maintain¬ 
ing u p /u* =2.00 and e = 0.95. Figure 6(a) | demonstrates 
that the distributions for density are the same after equal 
piston displacement. Scaling the mean kinetic energy 
(temperature) by the activation energy Ea = \u* 2 in 
Figure 6(a) we find the post-shock energy distributions 
are similar, tending towards a similar quasi-equilibrium 
state, where the kinetic energy tends to 5-8 % of the ac¬ 
tivation energy. This confirms that u p /u* is a scaling 
parameter for the dynamics. In our parametric study, 
we henceforth maintain u* = 10 and vary only u p and £. 
We also set the initial packing factor rji = 0.12, a param¬ 
eter we do not explore in the present study; see Sirmas 
et al. for its effect on the shock jump conditions in the 
case of non-dissipative collisions [53] , 


Dependence on e 

The role of £ on the shock structure is to control the 
relaxation rate. Figure [9(aJ| and [9(b)| show the distribu¬ 
tions of density and kinetic energy, respectively, for vary¬ 
ing £ and Up/u* = 2.00. Results show that decreasing e 
causes the kinetic energy to be excited and relaxed over a 
shorter length. This leads to a larger density gradient for 
lower £. The peak temperature decreases as £ decreases, 
owing to the increased dissipation during the initial ex¬ 
citation. The quasi-equilibrium states at the piston face 
show that the kinetic energies are similar, equal to ap¬ 
proximately 5% of the activation energy for £ = 0.80 and 
8% for £ = 0.95. This lower kinetic energy for decreasing 
£ leads to a somewhat higher density at the piston face 
after equal piston displacement. 

These trends are also seen by tracking the evolution of 
shock front, as shown in Figure [lO] Results show that 
decreasing £ generates a more rapid decay of the shock 
front. This is shown by the shock pulling towards the 
piston after a shorter time. These shocks are also closer 
to the piston, representing a more tightly packed relaxing 
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(a) (b) 

FIG. 6: Ensemble and coarse grained one-dimensional shock structure (a) of density and (b) kinetic energy after a 
piston displacement of x p = 138.7Ai for varying u p and u* with u p /u*= 2, and e = 0.95. 





FIG. 7: (Colour online) Comparison of the shock evolution for temperature (top) and pressure (bottom) obtained 
from MD for u p /u* = (a) 1.0, (b) 1.5, and (c) 2.0, with u* = 10, e = 0.95 and r)\ = 0.012. Selected particle paths 
(white) and forward running characteristics (black) are shown, where £ = x(t = 0). 


region. Although shocks develop faster with decreasing 
e, all shocks tend to approximately the same developed 
velocity. 

To conclude the parametric study, we look at the devel¬ 
oped shock morphology and variation of shock instability 
for varying u p /u* and e. These results are shown in Fig¬ 


ure 


11 


for u p /u* ranging from 1.00 to 3.00, and e of 0.80, 


0.90 and 0.95. The morphologies are taken after equal 
piston displacements of x p = 156.0Ai. Results show that 
the instabilities become prominent for all e with increas¬ 
ing u p /u*. As u p /u* increases, the frequency of these 
clusters extending from the equilibrium zone increases. 
























































(a) (b) 


FIG. 8: Ensemble and coarse grain averaged one-dimensional shock structure of (a) density and (b) temperature 
after a piston displacement of x p = 138.7Ai for different values of u p /u*, with £=0.95. 



(a) 
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FIG. 9: Ensemble and coarse grain averaged one-dimensional shock structure of (a) density and (b) temperature 
after a piston displacement of x p = 156.OAi for different values of e , with u p /u*= 2.0. 



FIG. 10: Evolution of shock front for varying e, with 
u* = 10 and u p = 20. 


The number of these instabilities also increases with de¬ 
creasing e. We find that the wavelength of these instabil¬ 
ities is on the same order as the relaxation length scales, 
as seen in distributions presented in Figs. [8] an m From 
these results, we see that the instabilities are noticeable 
for u p /u* > 1.00, with u p /u* = 1.00 difficult to discern, 
although this may be an artifact of the domain size. 


C. End States 

The variation of the end states for different shock 
strengths provides the shock Hugoniot, which can be 
used to assess whether the shock is unstable via the BZT 
and/or the DK instabilities, as discussed in the introduc¬ 
tion. Figure [12] shows the Hugoniot curve, on a pressure- 
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FIG. 11: Comparison of shock morphology for single realizations after a piston displacement of x p = 156.0Ai for 

different values of u p /u* and e, where u* = 10 and rji = 0.012. 
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FIG. 12: Results for Shock Hugoniot of equilibrium 
states, plotted with the elastic Hugoniot for 771 = 0.012, 
and the isotherm corresponding to u* = 10 . 


FIG. 13: Results for Rayleigh line (dashed) for 
Up/u* = 2.0 and £ = 0.95 plotted with the Shock 
Hugoniot (solid). 


specific volume (pv) plane, for the case studied of u* = 10 
and for £ = 0.95. Each point was evaluated in the post 
shock medium near the piston. As discussed above, the 
post shock state varies very slowly after the main re¬ 
laxation region (see for example Figure [ 8 (b)[ ), since an 
exponentially small fraction of the collisions remain in¬ 
elastic. For reference, we thus register the shock state as 
the point where the kinetic energy is 8 % of the activation 
energy. 

Results show that at sufficiently small piston speeds, 
i.e., Up/u* < 0.2, the post-shock state follows the theo¬ 
retical Hugoniot expected for a system of elastic disks, 
derived using Helfand’s equation of state [33] : 

P2 f( 1 -g) + ( 1 -’n > 2 (?) 

S( 1 -S".) 2 -|(i-s) 

where the jump in specific volume V 2 /V 1 = p\/P 2 - 

A transition occurs at approximately u p /u* = 0.2—0.3, 
corresponding to a high enough piston velocity activat¬ 
ing the inelastic collisions. Above this transition, the 
final state lies along the isotherm set by the activation 
threshold. Using Helfand’s equation of state [32 for the 
desired isotherm, here taken as uf. ms /u * 2 = 0.08, the 
final pressure is given by: 


P2 

Pi 


0.08— u* 2 
V2 



( 8 ) 


The evolution of the state from initial to final state 
across the steady shock is the so-called Rayleigh line. For 
further reference in our discussion of stability, Figure |T3| 
shows this path for the unstable case of u p /u* = 2.0 and 
£ = 0.95. 


3 

\ 

Q 


\ Isotherm velocity ratio 
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1 -- _ - - 
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u p /u* 


FIG. 14: Relationship between velocity of shock wave D 
and piston velocity u p , for u* = 10, £ = 0.90, and 
771 = 0.012. 


The speed of the shock waves were also determined by 
tracking the displacement of the shock front over sub¬ 
sequent time intervals. Figure [14] shows an example of 
the results for the shock velocity D for different values of 
u p /u* and £ = 0.90. Results show that at the lower veloc¬ 
ities, up to u p /u* = 0 . 2 , the velocities of the shock waves 
agree with the velocity predicted for elastic hard disks 
[33] ■ The shock velocity then deviates from this ideal 
behaviour between u p /u* = 0.3 — 1.0 until the velocity 
approaches D/u p « 1.0. The shock speed is in agree¬ 
ment with our theoretical prediction obtained by solving 
the jump equations for mass and momentum with the 
condition of isothermicity (Eq. Q). The shock velocity 
is well predicted by this solution for u p /u* > 0.3. 

To explain this transition occurring at u p /u* = 0.2 — 
0.3 we calculate the fraction of impact energy involved 
in the activated collisions, assuming a Boltzmann distri¬ 
bution for the state immediately behind the shock front. 
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This is completed by following the approach used in ki¬ 
netic theory to treat binary collisions, where one can be¬ 
gin with the rate of binary collisions per unit volume, 
written as [29j : 

6XP { ~^kT } C ° S ( 9 ) 

This term gives the rate of binary collisions of a system 
of disks of mass m with a number density n that have 
a relative speed in the range of g to g + dg 1 and an an¬ 
gle between the relative velocity and the line of action in 
the range of ip to ip + dip. The impact velocity, as men¬ 
tioned in Eq. <[T|) as the normal component of the relative 
velocity, is g n = geos ip. 

Multiplying Eq. § by ( g n ) 2 = (g cosip) 2 , and inte¬ 
grating over a range of g n , yields the energy along the 
line of action for collisions with impact velocities within 
this range of g n . Integrating g n from 0 to oo recovers the 
energy along the line of action for all collisions. Integrat¬ 
ing g n from u* to oo yields the energy seen along the line 
of action for impact velocities exceeding u*. From these 
results, we can calculate the fraction of the average en¬ 
ergy seen along the line of action for activated collisions, 
compared to that of all collisions. Acknowledging that 
u rms = 2fcT/ to, this ratio may be written as: 




(ff" 


= eX Pi -o 


g n > o 


2 v 2 
^ U'rms 


i + l u 


2 v 2 

^ a rms 


( 10 ) 


To evaluate the difference in this ratio for u p /u* = 0.2 
and 0.3, we assume that the temperature at the shock 
jump, before noticeable dissipation, can be estimated 
from elastic theory |33], where u 2 m . 


> u 2 . Using this 
equality in Eq. (flOl) allows us to approximate the frac¬ 


tion of impact energy that is sufficient to activate an in¬ 
elastic collision. The result for this ratio near the range 

As can be 


0 /u* = 0.2 and 0.3 is shown in Figure 


15 


seen, the fraction of impact energy that is activated is 
negligible for u p /u* = 0.2 (0.005 %) compared to that 
observed for u p /u* = 0.3 (2.5%). This clearly shows 
that Up/u* = 0.2 is not sufficiently strong to activate 
a significant number of inelastic collisions, and may be 
approximated using elastic jump conditions. However, 
u p /u* > 0.2 is shown to activate a more distinguishable 
number of collisions, which explains the transition from 
elastic theory seen around this value in the simulations. 


IV. DISCUSSION ON INSTABILITY 
MECHANISM 

A. Analysis of shock Hugoniot 

In the previous section, simulations showed that a 
shock structure does indeed become unstable with the 



u p /u* 


FIG. 15: Relationship between u p /u* and the fraction 
of the impact energy involved in activated collisions 
behind the shock front, assuming elastic jump 
conditions across shock wave. 


presence of dissipative collisions. Standard explanations 
for shock instability are related to the shock Hugoniot 
mmm- For the D’Yakov Kontorovich (DK) instabil¬ 
ity, the end states lying along sections of the Hugoniot 
having a positive slope are expected to have a corrugation 
type instability m■ Figure [12] shows that the Hugoniot 
does not take that form, ruling out the DK instability as 
an influencing mechanism. 

Another possible mechanism is if the fluid is of the 
BZT type or undergoes phase transitions. Shock split¬ 
ting is expected when the Rayleigh line, representing the 
state across the shock wave, intersects multiple points 
on the Hugoniot [34! ■ Such a behavior is possible near 
the transition u p /u* = 0.2 — 0.3 where the end state 
switches from lying on the elastic Hugoniot to lying on 
the isotherm. However, results demonstrate that it is for 
greater values of u p /u* that the shocks become unsta¬ 
ble. As seen in Figure 
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for Up/u* = 2.0, the Rayleigh 
line is far from this transition and does not intersect the 
Hugoniot in multiple locations, thus ruling out the insta¬ 
bility associated with shock splitting. Therefore, these 
mechanisms can be ruled out. 


B. Relaxation Rates and Comparison with 
Clustering Instability 

We now turn to another mechanism for instability pre¬ 
viously documented for homogeneous granular gases: the 
clustering instability in granular gases [23]. We wish to 
compare the residence time of the fluid in the shock struc¬ 
ture and the time scale required for clusters to develop 
within that element of fluid. Instability would ensue 
if the fluid resides within the relaxing region for longer 
times than required to develop the instability. 

The investigation of the clustering instability available 
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in the literature is for a homogeneous fluid at rest, which 
starts cooling while kept at constant volume. The evolu¬ 
tion of temperature before clustering is given by Haff’s 
law (see, for example, J3Q])- The parameters controlling 
this instability have been well documented [25H25] . and 
are not within the scope of the current work. One con¬ 
clusion we will adapt is that of Mitrano et al .: the onset 
of sensible clustering occurs when the evolution of gran¬ 
ular temperature deviates by 5% from Haff’s law j23] . 
Therefore by simulating the set of parameters observed 
in the shock waves, we can obtain the times scales for 
clustering necessary for comparison. 

To make a comparison between the instability of the 
constant specific volume case and the shock case, we com¬ 
pare the time evolution along the particle paths travers¬ 
ing the shock wave structure with the time history of 
cooling in a constant specific volume material element. 
For a meaningful comparison, this is done on time scales 
corresponding to the frequency of collisions, i.e., the lo¬ 
cal mean free time. This permits to automatically avoid 
accounting for density changes in calculating time scales. 
We adopt the same criterion for onset of instability as 
Mitrano et al ., and pose the question: How many local 
mean free times are required for the gas to develop in¬ 
stability ? and How many local mean free times does the 
shock transition last ? The comparison between these two 
time scales would permit to address whether the cluster¬ 
ing instability plays an important role. 

To obtain the characteristic time of clustering in terms 
of local mean free times, we first express Haff’s law in a 
time coordinate normalized by the local mean free time. 
Haff’s law expressed with time normalized by the initial 
mean free time, ti, may be expressed as [50] : 

m = _i_ (11) 

T ‘ (l+t|(l-e 2 )(l + 4«)) 2 


where 

16 (1 — e) (l — 2e 2 ) 
“ 2 ~ 57- 25e + 30e 2 (l-e) 


( 12 ) 


The relation between local and initial mean free time 
can be shown to be 


n _ Ai/u rm «(i) _ pb 2 (p) V-rms 

T A /u rms P\b 2 {r]i) ^rms(l) 


(13) 


Since the density p and packing factor p remain constant, 
(13) simplifies further to n/r = \/t]T\. 

Using this change in time scales in (11), we can obtain 
an expression for the theoretical evolution of temperature 
for a cooling homogeneous granular gas, in terms of time 
scaled by the local mean free time, i.e., t' = A. 

Constant volume clustering simulations were then con¬ 
ducted to determine the time when the energy of the 
system departs by more than 5% from Haff’s law, de¬ 
noting the time for the onset of clustering T c i ust . Since p 


varies across the shock structure, packing factors ranging 
from p = 0.05 — 0.25 were investigated using EDMD for 
e = 0.8, 0.9 and 0.95, and N =10,000. 

We now turn to establishing the relaxation time scale 
of the shocks. We track temperature along particle paths 
traversing the shock structure, with time integrated by 
using ( p~3| ) (select particle paths shown in Figure [3]). The 
relaxation time tr for each fluid element is obtained by 
fitting the temperature decay to an exponential decay 
equation: 


^=Hexp(--^ )+6 


(14) 


Figure 16 shows the results of tr/t for each particle 
element with varying u p /u* and e. The particles gener¬ 
ally experience fewer local mean free times to relax when 
e decreases or u p /u* increases. There are variations in 
relaxations times seen during the evolution of the shock 
wave. 

Since the specific particle paths along which the in¬ 
stability is triggered is unknown, we compute the mean 
value of tr/t for each set of parameters, as shown in Fig¬ 
ure [TT] The results show that at higher shock strengths 
(higher u p /u*) the time constant approaches some lim¬ 
iting value for each e. Given these results we now have 
a time scale to compare with the time to clustering in¬ 
stability T c i ust / t. Since the density increases across the 
shock wave, the value of p which contributes to the onset 
of instability can not be determined accurately. For this 
reason, the full range of clustering time for the range of 
p = 0.05 — 0.25 is compared. 

The results shown in Figure [17] indicate that there is 
practically no correlation between the observed shock in¬ 
stability and a residence time criterion. Unstable shocks 
are generally observed when the shock relaxation time is 
shorter than the clustering time. Likewise, stable shocks 
are observed when the shock relaxation time tr is longer 
than the characteristic time for clustering. There is al¬ 
most a perfect anti-correlation, suggesting that there is 
never sufficient time for a particle of fluid to develop a 
cluster, as it traverses the shock thickness during its re¬ 
laxation process. The results indicate that the clustering 
instability may not be the mechanism controlling shock 
instability. 


C. Role of initial transients on instability 

The instability of the shock was correlated above with 
the propensity of the relaxing medium to experience a re¬ 
pressurization event within the shock structure, as shown 
in Figure [7] For sufficiently small piston velocities (e.g. 
u p /u*= 1.0) the shock wave experiences a gradual decay 
in strength before attaining a developed structure propa¬ 
gating at a constant velocity. For this shock strength we 
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Up / u*=1 .0 + u p /u*=2.5 □ 

u D /u*=l.b x u o /u*=3.0 ■ 

Up/u*=2.0 x 

(a) € = 0.8 



% 


i/ p /i/*=1.0 + u p /u*=2.5 □ 

u D /u*= 1.5 x lL/u*= 3.0 ■ 

Up/u*= 2.0 x 

(b) e = 0.9 



u p /u'= 1.0 + u„/u*= 2.5 □ 

</,/u*=1.5 x </,/u*= 3.0 ■ 

Up/u‘= 2.0 x 

(c) e = 0.95 


FIG. 16: Exponential time constant of cooling experienced by shocked particle paths for different values of u p /u* 

and e, where £ = x(t = 0). 
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FIG. 17: Mean exponential time constant r# of shocked particle paths for different values of u p /u* and e, plotted 
with the range in time to clustering instability for similar values of e. Solid points represent simulations where 

unstable structures are seen. 


do not observe instabilities. However, as the piston ve¬ 
locity increases to u p /u*= 1.5 and above, the shock front 
stalls and pulls back towards the piston for a short period 
before attaining a developed structure. The evolution 
of these stronger shock waves exhibit a re-pressurization 
event experienced by the early particle paths. These pa¬ 
rameters also show the development of an unstable shock 
wave, suggesting a link between these initial transients 
and the stability. 

The re-pressurization during shock development sug¬ 
gests that the instability may be due to the pressure 
waves accelerating the flow along the piston. In this re¬ 
gion, very strong density gradients are established. These 
gradients become larger with increasing shock speed or 
decreasing s. These observations, and the type of in¬ 
stability observed with rolls forming along the density 
gradient, suggest that the mechanism controlling the in¬ 
stability is similar to Richtmyer-Meshkov or Rayleigh- 
Taylor type instabilities. It can be speculated that it is 
these wave phenomena that trigger multi-dimensional in¬ 
stabilities. This is also compatible with the absence of 


instability, other then the original pulsation, in ID simu¬ 
lations [7]. Further stability analysis of this initial tran¬ 
sient would be required, but its unsteadiness precludes 
using standard tools of linear analysis, such as the multi- 
mode approach. 

V. CONCLUSION 

The present study showed, for the first time, that re¬ 
laxing shock waves in granular gases develop instabilities, 
which take the form of convective rolls. Our investiga¬ 
tion of the possible mechanisms controlling the instabili¬ 
ties of shocks driven in relaxing media permitted to rule 
out several mechanisms. The reconstruction of the shock 
Hugoniot ruled out the D’Yakov-Kontorovich instability, 
as well as instability related to shock splitting. Results 
have shown shown that the shock waves develop the in¬ 
stability on similar times scales as the clustering instabil¬ 
ity seen in cooling granular gases. However, away from 
the stability limit, the time expected for clustering to oc- 
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cur is found to be always larger than the time scale for 
relaxation across the shock, suggesting that clustering 
instability is not the dominant mechanism. 

Nevertheless, the onset of instability was identified 
during the early stages of shock development and to 
correspond to the sufficient condition of an internal re¬ 
pressurization of the medium and subsequent pressure 
wave interaction with the density gradient. This sug¬ 
gests that the instability is of the Richtmyer-Meshkov 
type. Further study is required to quantify the interac¬ 
tions. 
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